/*
  Copyright 2010-2012 Patrik Jonsson, sunrise@familjenjonsson.org

  This file is part of Sunrise.

  Sunrise is free software; you can redistribute it and/or modify
  it under the terms of the GNU General Public License as published by
  the Free Software Foundation; either version 3 of the License, or
  (at your option) any later version.

  Sunrise is distributed in the hope that it will be useful,
  but WITHOUT ANY WARRANTY; without even the implied warranty of
  MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
  GNU General Public License for more details.

  You should have received a copy of the GNU General Public License
  along with Sunrise.  If not, see <http://www.gnu.org/licenses/>.

*/

/// \file
/// Contains classes related to emission of rays from Arepo cells.

#ifndef __arepo_emission__
#define __arepo_emission__

#include "emission.h"
#include "arepo_grid.h"
#include "almost_equal.h"

namespace mcrx {
  template <typename, typename> class arepo_cell_emission;
}

#ifdef WITH_AREPO

/** Represents volumetric emission from an Arepo Voronoi cell. The
    position within the Voronoi cell is drawn randomly from a
    distribution proportional to the density within the cell using
    rejection sampling. The set up of the bounding box information is
    a bit problematic as this requires the circumsphere centers of the
    Delaunay tessellation, which is not stored indexed by the
    cells. For this reason, after creating all the cell_emission
    objects, a final pass is done which loops over the Delaunay
    tetrahedrons and updates the bounding box information for the
    cells associated with that tetrahedron. This is done by the static
    post_setup() function, which must be called after the grid has
    been constructed. This is unfortunate, but the alternative is to
    loop over the DT array once per created cell, which would be
    prohibitively expensive.  */
template <typename chromatic_policy, 
	  typename rng_policy_type>
class mcrx::arepo_cell_emission :
  public volumetric_emission<chromatic_policy, rng_policy_type> {
public:
  typedef volumetric_emission<chromatic_policy, rng_policy_type> T_base;
  typedef typename T_base::T_lambda T_lambda;
  typedef typename T_base::T_biaser T_biaser;
  typedef typename T_base::T_ray T_ray;
  typedef rng_policy_type T_rng_policy;

  /** To avoid circularly dependent template parameters, it is
  hardcoded that the grid cells contain this class.  This precludes
  making grids which contain objects which in turn contains
  grid_cell_emission objects, but was chosen for simplicity and the
  lack of the current need for anything more advanced. */
  typedef arepo_grid_cell<cell_data<arepo_cell_emission, absorber<T_lambda> > > T_cell;

private:
  /** True if the bounding box setup has been done. This is just used
      to check that the post_setup step has been done. */
  bool setup_done_; 

  const T_cell* cell_; ///< Pointer to the cell doing the emission.

  T_float rho_min_; ///< The mininimum density in the cell. In Arepo units.
  T_float rho_max_; ///< The maximum density in the cell. In Arepo units.

  vec3d min_; ///< The min coordinates of the cell bounding box. In Arepo units.
  vec3d max_; ///< The max coordinates of the cell bounding box. In Arepo units.

public:
  /** Default constructor creates an empty/invalid emitter. */
  arepo_cell_emission() : 
    volumetric_emission<chromatic_policy, rng_policy_type>(), 
    setup_done_(false), cell_(0) {};

  /** Constructor for polychromatic situations, which needs to know
      the wavelength vector. */
  template<typename T_grid_iterator>
  arepo_cell_emission (const T_lambda& em, const T_lambda* lptr,
		       const T_grid_iterator& c):
    volumetric_emission<chromatic_policy, rng_policy_type>(em, lptr), 
    setup_done_(false), cell_(&(*c)),
    rho_min_(SphP[cell_->index()].Density*arepo::mcon/(arepo::lcon*arepo::lcon*arepo::lcon)),
    rho_max_(rho_min_),
    min_(blitz::huge(T_float())), max_(blitz::neghuge(T_float()))
  {
    // Note that the bounding box can not be set correctly at this
    // point. These must be updated by looping over the tetrahedrons
    // that are connected to this mesh point.

    // we are initialized with invalid emission
    //ASSERT_ALL(em==em);
    //ASSERT_ALL(em<1e300);
  };

  /** This constructor is for non-polychromatic situations where we
      don't have a concept of wavelength. */
  template <typename T_grid_iterator>
  arepo_cell_emission (const T_lambda& em, const T_grid_iterator& c):
    volumetric_emission<chromatic_policy, rng_policy_type>(em), 
    setup_done_(false), cell_(&(*c)),
    rho_min_(SphP[cell_->index()].Density*arepo::mcon/(arepo::lcon*arepo::lcon*arepo::lcon)),
    rho_max_(rho_min_),
    min_(blitz::huge(T_float())), max_(blitz::neghuge(T_float()))
  {
    // Note that the bounding box can not be set correctly at this
    // point. These must be updated by looping over the tetrahedrons
    // that are connected to this mesh point.

    // we are initialized with invalid emission
    //ASSERT_ALL(em==em);
    //ASSERT_ALL(em<1e300);
  };

  // copy constructor and assignment operator are implicitly generated

  /** This accessor is so that a grid can use emitters both directly
      and with the the cell_data class, while still retaining that
      possibility. Hopefully use of absorber gets optimized out...  */
  arepo_cell_emission& get_emitter() {return *this;};
  const arepo_cell_emission& get_emitter() const {return *this;};

  virtual ray_distribution_pair<T_lambda> emit (T_biaser, T_rng_policy&, 
						T_float&) const; 

  boost::shared_ptr<emission<chromatic_policy, rng_policy_type> > clone() const {
    return boost::shared_ptr<emission<chromatic_policy, rng_policy_type> >(new arepo_cell_emission<chromatic_policy, rng_policy_type>(*this));};

  /** Update the bounding box information with the fact that p is a
      point on the cell face. P should be in Arepo units. */
  void update_bounding_box(const vec3d& p) {
    min_=where(p<min_, p, min_); max_=where(p>max_, p, max_);
    // update density min and max values with the values at this location
    const T_float rho = arepo::cell_density(cell_->index(), p);
    assert(rho>=0);
    rho_min_=std::min(rho_min_, rho);
    rho_max_=std::max(rho_max_, rho);
 };

  /** This function sets the bounding boxes of the cells by looping
      over all the Arepo tetrahedrons and updating the bounding boxes
      with the positions of the circumsphere centers of tetrahedrons
      connected to the mesh point. */
  template<typename T_data>
  static void update_bounding_boxes(const arepo_grid<T_data>& g);

  /** This routine uses the arepo_grid object to setup of the bounding
      box information. It MUST be called before emitting rays. */
  template<typename T_data>
  static void post_setup (const arepo_grid<T_data>& g) {
    update_bounding_boxes(g); };
};


/** Emit a ray from a location within the cell drawn based on the
    density distribution. The direction is drawn isotropically. */
template <typename chromatic_policy, 
	  typename rng_policy_type>
mcrx::ray_distribution_pair<typename mcrx::arepo_cell_emission<chromatic_policy, rng_policy_type>::T_lambda>
mcrx::arepo_cell_emission<chromatic_policy, rng_policy_type>::
emit (T_biaser, T_rng_policy& rng, T_float& prob) const
{
  using namespace arepo;
  assert(setup_done_);

  // we draw a point consistent with the density distribution in the
  // cell using rejection sampling
  const int sph_idx = cell_->index();

  vec3d pos;
  int n=0;

  DEBUG(1,std::cout << "arepo_cell_emission: Drawing emission point for cell " << sph_idx << endl;);

  while(true) {
    ++n;
    DEBUG(1,std::cout << "\tAttempt " << n << "\n";);
    if((n>100) && (n%100==0)) {
      std::cout << "Warning: many rejections in arepo_emission:" << n
		<< std::endl;
    }

    assert(n<10000);

    // first we draw a point from within the bounding box (which is in
    // Arepo units)
    blitz::TinyVector<T_float, 3> r;
    rng.rnd(r);
  
    pos = (max_-min_)*r + min_;

    if(!point_inside(sph_idx, pos)) {
      // If it's outside we try a new position
      DEBUG(1,std::cout << "\tPoint " << pos << " outside cell, rejecting\n";);
      continue;
    }

    DEBUG(1,std::cout << "\tPoint " << pos << " inside cell, proceeding...\n";);

    // if we are inside, we now accept with a probability of rho(pos)/rho_max.
    const T_float r2(rng.rnd());
    if (r2*rho_max_ < rho_min_) {
      // we can accept without testing, because we know the density in
      // this spot is larger than rho_min_.
      DEBUG(1,std::cout << "\taccepting density without testing\n";);
      break;
    }

    // we could not accept without testing, so we need to calculate
    // the actual density at this point
    const T_float rho = arepo::cell_density(sph_idx, pos);
    if (r2*rho_max_ < rho) {
      // we accept this point
      DEBUG(1,std::cout << "\tDensity fraction " << r2 << " accepted.\n";);
      break;
    }

    DEBUG(1,std::cout << "\tDrawn density fraction " << r2 << " higher than " << rho/rho_max_ << ", rejecting\n";);
  }

  // we come out here with a point correctly drawn according to the
  // density distribution in the cell. Convert the point to a Sunrise
  // position.
  pos = periodic_wrap_point(pos)*lcon;

  // Now draw an angle and we're done.
  blitz::TinyVector<T_float, 2> r;
  rng.rnd(r);

  // and the emission direction (from an isotropic distribution)
  const T_float dphi = r [0]*2*constants::pi; 
  const T_float dcos_theta = 2*r [1] - 1;
  const T_float dsin_theta = sqrt (1- dcos_theta*dcos_theta);

  // probability is density/mass from position, and 1/4pi from
  // direction. calculate if we ever need it.
  prob = blitz::quiet_NaN(T_float());

  //and create the ray
  return ray_distribution_pair<T_lambda> 
    (boost::shared_ptr<T_ray> 
     (new T_ray (pos,
		 vec3d (dsin_theta*cos (dphi), dsin_theta*sin (dphi), 
			dcos_theta),
		 initialize(this->emission_, 1.))), 
     this->emission_,
     this->ad_);
}

template <typename chromatic_policy, 
	  typename rng_policy_type>
template <typename T_data>
void
mcrx::arepo_cell_emission<chromatic_policy, rng_policy_type>::
update_bounding_boxes(const arepo_grid<T_data>& g)
{
  using namespace arepo;
  tetra* DT = T->DT;
  point* DP = T->DP;
  tetra_center* DTC = T->DTC;
  std::cout << "Calculating bounding boxes of cells." << std::endl;

  // Start by looping over all tetras
  for(int t=0; t<T->Ndt; ++t) {
    if( DT[t].t[0]<0 || 
	DT[t].p[0] == DPinfinity || DT[t].p[1] == DPinfinity || 
	DT[t].p[2] == DPinfinity || DT[t].p[3] == DPinfinity )
      // tetras with t[0]<0 or that contain an infinity point do not
      // have computed circumspheres. (It's fine, <0 tetras are
      // deleted and the infinity points are well outside our box.)
      continue;

    const vec3d dtc (vec3d(DTC[t].cx, DTC[t].cy, DTC[t].cz));

    // For each tetrahedron, we loop over the four vertices
    for(int p=0; p<4; ++p) {

      // find the voronoi cell connected to this vertex. Here, we need
      // only concern ourselves with primary mesh points. They may be
      // duplicated as ghosts, but to get all dtc's of a primary mesh
      // point, we need only process tetras that are connected to the
      // primary point itself.  to the primary mesh point
      const int dp_idx =DT[t].p[p];

      if( (DP[dp_idx].task==ThisTask) && 
	  (dp_idx>=0) && (dp_idx<N_gas) &&
	  (SphP[dp_idx].Density>0) ) {

	arepo_cell_emission& c = g.cell(dp_idx).data()->get_emitter();
	// and update bounding box, taking care to wrap the point
	// around the boundary
	const vec3d center(SphP[dp_idx].Center[0],
			   SphP[dp_idx].Center[1],
			   SphP[dp_idx].Center[2]);
	c.update_bounding_box(arepo::periodic_wrap_point(dtc, center));
	c.setup_done_=true;
      }
    }
  }
}


#endif
#endif
